#!/usr/bin/perl -w
=head USAGE

  panda_enter_PFAM_HMM_into_postgres.pl

=head PURPOSE

    This script takes
       1) an HMM file containing multiple HMMs
       2) the PFAM data file (Pfam-A.seed)
    and enters them into tables of a specified postgres database

=head REMARKS

=cut
#die and warn check next

use strict;
use Getopt::Long;
use warnings;
use POSIX qw(:float_h); #for DBL_MAX
use DBI;
use DBD::Pg;
use String::CRC;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

my $usage = <<'USAGE';

USAGE:

    panda_enter_PFAM_HMM_into_postgres.pl
                            --align_data <pfam_file>
                            --hmm_file   <pfam_hmm>
                            [--dbname   <panda>]
                            [--user      <any_db_username>]
                            [--password  <secret>]
                            [--help]

  				
  				
    This script takes
       1) an HMM file <pfam_hmm> containing multiple HMMs
       2) the PFAM data file <pfam_file> (e.g. Pfam-A.seed)
    and enters them into tables of a specified postgres database

    The PFAM seed alignment file <pfam_file> provides
    the accession, name, description, author and
	alignment method for each PFAM domain
	
	
USAGE




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# options variables
my $help;
my $dbname		='panda';
my $username	='';
my $userpassword='';
my $pfam_file	='';
my $hmm_file	='';

GetOptions (
			'help'			=> \$help,

			'dbname=s'		=> \$dbname,
			'user=s'		=> \$username,
			'password=s'	=> \$userpassword,

			'hmm_file=s'	=> \$hmm_file,
			'align_data=s'	=> \$pfam_file,
			);



die $usage if ($help);
die "\nERROR:\n\tMissing Argument: (align_data)\n" .$usage
				unless ($pfam_file);
die "\nERROR:\n\tMissing Argument: (hmm_file)\n" .$usage
				unless ($hmm_file);






#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Various functions

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

sub create_postgres_tables($);
sub parse_pfam_domains_data($\%);
sub parse_and_save_HMM_file($$\%\%\@);
	sub lookup_domain_data_for_HMM(\%\%);
		sub bitscore_to_P_value($$$);
	sub send_to_db($$\%);
		sub remove_matching_panda_domain_data($$);
sub delete_obsolete_HMMs(\%$);
sub drop_indices($$);




#_________________________________________________________________________________________

#	drop_indices($$)

#		drops all indices for a given table and
#		returns the commands needed to recreate them

#_________________________________________________________________________________________
sub drop_indices($$)
{
	my ($dbh, $table_name) = @_;

	# get all index names and the commands to recreate them
	my $sth_get_indices	= $dbh->prepare("SELECT indexname, indexdef FROM pg_indexes "
							."WHERE tablename = ? ");
	$sth_get_indices->execute($table_name) or die "\nError:\n$DBI::errstr\n";
	my @recreate_commands;
	my @index_names;
	while (my @res = $sth_get_indices->fetchrow_array())
	{
		push(@recreate_commands, $res[1]);
		push(@index_names, $res[0]);
	}
	
	# return nothing if no indices on this table
	return () if (@index_names == 0);
	
	# drop all indices
	for (@index_names)
	{
		$dbh->do("DROP INDEX $_")
				  or die "\nError:\n$DBI::errstr\n";
	}
	
	return @recreate_commands;
}



#_________________________________________________________________________________________


#	create_postgres_tables


#_________________________________________________________________________________________

sub create_postgres_tables($)
{
	my ($dbh) = @_;	
	# this procedure tends to fail (if the tables already exist)
	# so no error messages but auto commit
	my @original_dbh_flags = @{$dbh}{'PrintError', 'AutoCommit', 'RaiseError'};
	@{$dbh}{'PrintError', 'AutoCommit', 'RaiseError'} = (0, 1, 0);
	if ($dbh->do( 'CREATE TABLE domains ('.
					'hmm_name          TEXT,'.				#	NAME
					'dom_code          INTEGER,'.			#	ACC
					'thres_true        DOUBLE PRECISION,'.	#	TC
					'thres_maybe       DOUBLE PRECISION,'.	#	NC
					'thres_repeat      DOUBLE PRECISION,'.	#	THRES_REPT
					'thres_subfamily   DOUBLE PRECISION,'.	#	THRES_FAM
					'definition        TEXT,'.				#	DEFN
					'description       TEXT,'.				#	DESC
					'binding_catalysis TEXT)'			)	#	BIND_CATA
		)
	{
		$dbh->do( 'CREATE INDEX i_domains__dom_code on domains(dom_code)');
	}

	$dbh->do( 'CREATE TABLE archived_domains ('.
					'hmm_name          TEXT,'.				#	NAME
					'dom_code          INTEGER,'.			#	ACC
					'thres_true        DOUBLE PRECISION,'.	#	TC
					'thres_maybe       DOUBLE PRECISION,'.	#	NC
					'thres_repeat      DOUBLE PRECISION,'.	#	THRES_REPT
					'thres_subfamily   DOUBLE PRECISION,'.	#	THRES_FAM
					'definition        TEXT,'.				#	DEFN
					'description       TEXT,'.				#	DESC
					'binding_catalysis TEXT,'.				#	BIND_CATA
					'dom_date_archived TIMESTAMP '.
							'DEFAULT CURRENT_TIMESTAMP)'	);

	if ($dbh->do('CREATE TABLE hmms ('.
					'hmm_name          TEXT,'.				#	NAME
					'hmm_acc           TEXT,'.				#	ACC
					'dom_code          INTEGER,'.			#	ACC
					'hmm_file          TEXT,'.				#	FILENAME
					'hmm_crc1          OID,'.				#	CRC1
					'hmm_crc2          OID,'.				#	CRC2
					'hmm_db_origin     TEXT,'.				#	'PFAM'
					'aln_method        TEXT,'.				#	ALN_METHOD
					'db_state          TEXT,'.				#	DB_STATE
					'author            TEXT,'.				#	AUTH
					'hmm               TEXT,'.				#	HMM			
					'hmm_date_added    TIMESTAMP '.
							'DEFAULT CURRENT_TIMESTAMP)'	)
		)
	{
		$dbh->do( 'CREATE INDEX i_hmms__dom_code on hmms(dom_code)');
	}
	$dbh->do('CREATE TABLE new_hmms ('.
					'dom_code          INTEGER)'				);
	$dbh->do('CREATE TABLE archived_hmms ('.
					'hmm_name          TEXT,'.				#	NAME
					'hmm_acc           TEXT,'.				#	ACC
					'dom_code          INTEGER,'.			#	ACC
					'hmm_file          TEXT,'.				#	FILENAME
					'hmm_crc1          OID,'.				#	CRC1
					'hmm_crc2          OID,'.				#	CRC2
					'hmm_db_origin     TEXT,'.				#	'PFAM'
					'aln_method        TEXT,'.				#	ALN_METHOD
					'db_state          TEXT,'.				#	DB_STATE
					'author            TEXT,'.				#	AUTH
					'hmm               TEXT,'.				#	HMM			
					'hmm_date_added    TIMESTAMP,'.
					'hmm_date_archived TIMESTAMP '.
							'DEFAULT CURRENT_TIMESTAMP)'	);
	if 	($dbh->do( 'CREATE TABLE domains_per_protein ('.
					'pi1               OID,'.
					'pi2               OID,'.
					'dom_code          INTEGER,'.
					'dom_name          TEXT,'.
					'start             INTEGER,'.
					'finish            INTEGER,'.
					'p_value           DOUBLE PRECISION,'.
					'confidence        INTEGER)'))
	{
		$dbh->do( 'CREATE INDEX i_domprot__pi on domains_per_protein(pi1,pi2)');
		$dbh->do( 'CREATE INDEX i_domprot__dom_code on domains_per_protein(dom_code)');
		$dbh->do( 'CREATE INDEX i_domprot__dom_name on domains_per_protein(dom_name)');
	}
	@{$dbh}{'PrintError', 'AutoCommit', 'RaiseError'} = @original_dbh_flags;
}	


#_________________________________________________________________________________________

#	parse_pfam_domains_data

#		Get pfam data
#		Parse the Pfam-A.seed file to get domain annotation which is not
#		in the HMM file itself:
	
#			i.e. description, author and alignment method

#		takes
#				1) the pfam file name and
#				2) a hash keyed on acc_name for returning the data

#_________________________________________________________________________________________
sub parse_pfam_domains_data($\%)
{
	my ($pfam_file, $pfam_domains_data) = @_;	

	# retrieve five data types from pfam data file
	my %PFAM_DATA_TYPES;
	@PFAM_DATA_TYPES{qw(ID AC DE AU AL)} = ();


	open PFAM_FILE, $pfam_file or die "Error:\nCould not open $pfam_file.\n$!\n";
	print STDERR '[';
	my %data;
	my $rec_cnt;
	while (<PFAM_FILE>)
	{
		# save data at the end of each record
		if (m{^//})
		{
			print STDERR "." if ($rec_cnt++ % 20 == 0);
			
			# check that all five data categories were saved
			die "$!\n" unless (keys %data) == 5;

			# save data
			$pfam_domains_data->{$data{'AC'}.'|'.$data{'ID'}} = [@data{qw(DE AU AL)}];

			# clear variables
			%data = ();
		}

		next unless /#=GF +([IAD][DCEUL]) +(.+)\n/;

		# check that this is the data type we want to save
		next unless exists $PFAM_DATA_TYPES{$1};

		$data{$1} = $2;
	}
	close PFAM_FILE;
	print STDERR "]\n$rec_cnt\t Pfam domain records parsed.\n";
}	


#_________________________________________________________________________________________

#	parse_and_save_HMM_file


#_________________________________________________________________________________________
sub parse_and_save_HMM_file($$\%\%\@)
{
	my ($hmm_file, $dbh, $HMMs, $pfam_domains_data, $obsolete_hmms) = @_;


	open HMM_FILE, $hmm_file or die "Error:\nCould not open $hmm_file.\n$!\n";

	# create statement handle for checking whether the entry already exists
	my $sth_chk_exists	= $dbh->prepare(
									'SELECT hmm_name, hmm_acc, hmm_crc1, '.
									'hmm_crc2 '.
									'FROM hmms '.
									'WHERE dom_code = ?');

	# fields to be parsed from the HMM file
	my %HMM_FILE_FIELDS;
	@HMM_FILE_FIELDS{qw(NAME ACC DESC TC NC EVD)} = ();




	# domain data parsed from the HMM
	my %HMM_fields;
	my $rec_cnt = 0;
	my $hmm_for_crc = '';

	print STDERR '[';
	while (<HMM_FILE>)
	{
		# concatenate separate lines of HMM together
		$HMM_fields{'HMM'} .= $_;
		$hmm_for_crc .= $_ unless (/^COM/ or /^DATE/);
		
		# end of record, process
		if (/^\/\//)
		{
			$HMM_fields{'FILENAME'}	= $hmm_file;
			# remove path
			$HMM_fields{'FILENAME'} =~ s[.*/][];

			my ($crc1, $crc2) = crc ($hmm_for_crc, 64);
			$HMM_fields{'CRC1'}	= $crc1;
			$HMM_fields{'CRC2'}	= $crc2;


			# make sure all fields have been parsed
			foreach (keys %HMM_FILE_FIELDS)
			{
				die "Error:\nFailed to parse a field [$_] from the HMM file $hmm_file\n"
												unless defined ($HMM_fields{$_});
			}
		
			# only save if can find pfam data corresponding to the HMM
			if (lookup_domain_data_for_HMM(%HMM_fields, %$pfam_domains_data))
			{
				$HMMs->{$HMM_fields{'DOM_CODE'}} = 0;
			

				# save hmm in list
				if (send_to_db($dbh, $sth_chk_exists, %HMM_fields))
				{
					$HMMs->{$HMM_fields{'DOM_CODE'}} = 1;
					print STDERR ":" if ($rec_cnt++ % 20 == 0);
				}
			}
			else
			{
				push(@$obsolete_hmms, $HMM_fields{'NAME'}.'['.$HMM_fields{'ACC'}.']');
				print STDERR "!" if (@$obsolete_hmms % 20 == 0);
			}

			# clear variables
			%HMM_fields = ();
			$hmm_for_crc = '';
			
		}
		# only save matching NAME ACC DESC TC NC EVD
		next unless /^([A-Z]+) +(.*)\n/;
		next unless exists $HMM_FILE_FIELDS{$1};
		$HMM_fields{$1} = $2;

	}
	close HMM_FILE;
	print STDERR "]\n$rec_cnt\t Pfam HMM records parsed.\n";
}	


#_________________________________________________________________________________________

#	lookup_domain_data_for_HMM

#				Takes data from an HMM record, prepares it for sending to the db

#		Parameters:



#_________________________________________________________________________________________
sub lookup_domain_data_for_HMM(\%\%)
{
	my ($HMM_fields, $pfam_domains_data) = @_;

	# use pfam data
	my $key = $HMM_fields->{'ACC'}.'|'.$HMM_fields->{'NAME'};
	return undef unless (exists $pfam_domains_data->{$key});
	
	# note that the pfam decription is the panda definition
	@{$HMM_fields}{qw (DEFN AUTH ALN_METHOD) } = @{$pfam_domains_data->{$key}};



	# calculate E values cut-offs from the bitscores
	my ($mu, $lambda) = split (/ +/, $HMM_fields->{'EVD'});
	delete $HMM_fields->{'EVD'};
	no warnings ('numeric');
	$HMM_fields->{'TC'} = bitscore_to_P_value(	$HMM_fields->{'TC'}+0,
												$mu,
												$lambda);
	$HMM_fields->{'NC'} = bitscore_to_P_value(	$HMM_fields->{'NC'}+0,
												$mu,
												$lambda);
	use warnings ('numeric');
	
	$HMM_fields->{'DB_ORIGIN'}= 'PFAM';
	
	# set fields absent from pfam to undef
	my @DOMAIN_FIELDS_NOT_IN_PFAM = qw (THRES_REPT THRES_FAM DESC DB_STATE BIND_CATA);
	@{$HMM_fields}{@DOMAIN_FIELDS_NOT_IN_PFAM} = ();

	# set domain code
	# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	# PFAM accessions start at 50000
	$HMM_fields->{'ACC'} =~ m/PF(\d+)/;
	$HMM_fields->{'DOM_CODE'} = $1 + 50000;

}	


#_________________________________________________________________________________________

#	bitscore_to_E_value

#				Calculates the E-value from the bitscore given an EVD fit

#		Parameters:

#				bitscore, mu, lambda, database_size


#_________________________________________________________________________________________
sub bitscore_to_P_value($$$)
{

	my($bit_score, $mu, $lambda) = @_;


	my $pval = 1.0 / (1.0 + (exp($bit_score * 0.69314718)));

	# use approximation if no EVD
	return ($pval)
			if(! defined($mu));

	#Rounding to  1.0
	return 1.0
			if(( $lambda * ($bit_score - $mu)) < -1.0 * log(log(DBL_MAX)));

	#Rounding to  0.0
	return 0.0
			if(( $lambda * ($bit_score - $mu)) > log(DBL_MAX));

	my $y = exp(-1.0 * $lambda * ($bit_score - $mu));
	if($y >= 1.0e-7)
	{
		$y = (1.0 - exp (-1.0 * $y));
	}

	if($y < $pval){ $pval = $y }

	return ($pval);

}

#_________________________________________________________________________________________

#	send_to_db

#				Enters domain data into the postgres table

#		Parameters:

#				DBI database handle
#				DBI handle for SQL statement checking whether the entry already exists
#				List of data in a hash
#				List of data keys for the hash


#_________________________________________________________________________________________

sub send_to_db($$\%)
{
	# get parameters
	(my ($dbh, $sth, $HMM_fields)) = @_;
	
	

	# make sure all fields to be saved in the HMM file have been parsed
	my @hmms_table = qw (NAME ACC DOM_CODE FILENAME CRC1 CRC2
							DB_ORIGIN ALN_METHOD DB_STATE AUTH HMM);
	my @domains_table = qw (NAME DOM_CODE TC NC THRES_REPT
							THRES_FAM DEFN DESC BIND_CATA);



	# make sure all fields are parsed
	foreach (@domains_table, @hmms_table)
	{
		die "$!\n" unless exists $HMM_fields->{$_};
	}




	# return if the HMM is already up-to-date in the database
	$sth->execute($HMM_fields->{'DOM_CODE'}) or die "$!\n" ;
	if ($sth->rows == 1)
	{
		my ($name, $acc, $crc1, $crc2) = $sth->fetchrow_array();
		$sth->finish();
		if ($acc	eq $HMM_fields->{'ACC'}		&&
			$name	eq $HMM_fields->{'NAME'}	&&
			$crc1	eq $HMM_fields->{'CRC1'}	&&
			$crc2	eq $HMM_fields->{'CRC2'})
		{
			return undef;
		}
	}
	else
	{
		$sth->finish();
	}
	
	
	# remove all domain data for this code
	remove_matching_panda_domain_data($dbh, $HMM_fields->{'DOM_CODE'});


	
	# Enter in new HMM data
	$dbh->do( 	'INSERT INTO domains VALUES ('.
				join (	', ', ('?') x @domains_table) . ')',
				undef, @{$HMM_fields}{@domains_table});
	# add quotes for text
	for (qw (NAME ACC FILENAME DB_ORIGIN ALN_METHOD DB_STATE AUTH HMM))
	{
		# get rid of nulls
		$HMM_fields->{$_} |= '';
		$HMM_fields->{$_} = $dbh->quote($HMM_fields->{$_});
	}
	for (qw (CRC1 CRC2))
	{
		$HMM_fields->{$_} = "OID(".$HMM_fields->{$_}.")";
	}
	my $sql = 'INSERT INTO hmms VALUES ('.
				join (	', ', @{$HMM_fields}{@hmms_table}).
				')';
	eval{
	$dbh->do($sql);
	};
	if ($@)
	{
		print $sql;
		exit;
	}


	# tag HMM as new entry
	$dbh->do( 	'INSERT into new_hmms VALUES (?)',
				undef, @$HMM_fields{'DOM_CODE'});

	return 1;
}				


#_________________________________________________________________________________________

#	remove_matching_panda_domain_data

#				Removes domain data in panda which matches a domain code

#		Parameters:

#				DBI database handle
#				domain_code


#_________________________________________________________________________________________
sub remove_matching_panda_domain_data($$)
{
	(my ($dbh, $dom_code)) = @_;

	# archive matching old HMMs
	$dbh->do(	'INSERT INTO archived_hmms ('.
				'SELECT * from hmms WHERE dom_code = ?)',
				undef, ($dom_code));
	$dbh->do(	'INSERT INTO archived_domains ('.
				'SELECT * from domains WHERE dom_code = ?)',
				undef, ($dom_code));

	# clear matching old HMMs
	$dbh->do( 	'DELETE FROM hmms WHERE dom_code = ?',
				undef, ($dom_code));
	$dbh->do( 	'DELETE FROM domains WHERE dom_code = ?',
				undef, ($dom_code));
	$dbh->do( 	'DELETE FROM new_hmms WHERE dom_code = ?',
				undef, ($dom_code));
			

	# clear matching HMM Seach results
	$dbh->do( 	'DELETE FROM domains_per_protein WHERE dom_code = ?',
				undef, ($dom_code));
}	
#_________________________________________________________________________________________


#	delete_obsolete_HMMs


# 		Delete obsolete PFAM HMMs
#		N.B. Calls bommit


#_________________________________________________________________________________________
sub delete_obsolete_HMMs(\%$)
{
	my ($HMMs, $dbh) = @_;
	my @curr_HMMs;

	# create statement handle for checking whether the entry already exists
	my $sth_get_all_domains = $dbh->prepare('SELECT dom_code FROM hmms WHERE '.
											"hmm_db_origin = 'PFAM'");
	$sth_get_all_domains->execute() or die "$!\n" ;
	while (my ($dom_code) = $sth_get_all_domains->fetchrow_array())
	{
		push (@curr_HMMs, $dom_code);
	}
	$sth_get_all_domains->finish();

	my @obsolete_HMMs;
	foreach(@curr_HMMs)
	{
		push(@obsolete_HMMs, $_) unless exists ($HMMs->{$_});
	}
	
	foreach my $dom_code(@obsolete_HMMs)
	{
		remove_matching_panda_domain_data($dbh, $dom_code);
	}

	
}	


















#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Core logic!!

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# connect to database
my $dbh = DBI->connect( "dbi:Pg:dbname=$dbname",
                        $username,
                        $userpassword,
						{
                          RaiseError => 1,
                          AutoCommit => 0,
						  PrintError => 1
                        }
                      )
		or die "Database connection not made: $DBI::errstr";

# make sure HMM / domain tables are present in the database
create_postgres_tables($dbh);

#drop all indices on the tables
my @sql_recreate_indices_new_hmms	= drop_indices($dbh, 'new_hmms');
my @sql_recreate_indices_domains	= drop_indices($dbh, 'domains');
my @sql_recreate_indices_hmms		= drop_indices($dbh, 'hmms');



# get the domain data for pfam
my %pfam_domains_data;
parse_pfam_domains_data($pfam_file, %pfam_domains_data);



		
# save list of PFAM HMMs
my %HMMs;

# parse the HMM file
#	calls send_to_db iteratively
my @obsolete_hmms;
parse_and_save_HMM_file($hmm_file, $dbh, %HMMs, %pfam_domains_data, @obsolete_hmms);
$dbh->commit();

				
# remove data for obsolete domains
delete_obsolete_HMMs(%HMMs, $dbh);
$dbh->commit();


# make sure new_hmm tables does not contain duplicates
$dbh->do("CREATE TABLE temp_new_hmms as SELECT DISTINCT on (dom_code) * from new_hmms");
$dbh->do("DROP TABLE new_hmms");
$dbh->do("ALTER TABLE temp_new_hmms rename to new_hmms");


#recreate indices
$dbh->do($_) or die "\nError:\n$DBI::errstr\n" for (@sql_recreate_indices_new_hmms);
$dbh->do($_) or die "\nError:\n$DBI::errstr\n" for (@sql_recreate_indices_domains);
$dbh->do($_) or die "\nError:\n$DBI::errstr\n" for (@sql_recreate_indices_hmms);


# commit changes and disconnect
$dbh->commit();
$dbh->disconnect();

# report results
my $cnt_HMMs;
$cnt_HMMs += $HMMs{$_} for(keys %HMMs);
print STDERR "\nSuccess!!!\n$cnt_HMMs PFAM HMMs and their data entered into panda.\n";

# list of hmms which are obsolete
if (@obsolete_hmms)
{
	print STDERR scalar @obsolete_hmms, " HMM records do not have ".
			"corresponding data in the PFAM domain annotations file.\n\n";
	print "These PFAM domains are most probably obsolete and will not be entered into panda:\n";
	print join "\n", @obsolete_hmms, "";
}
